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INTRODUCTION 

The term plasma-sheath problem denotes any boundary-value problem described by the 
Poisson and steady-state Boltzmann equations. The charge density is obtained through 
solutions of the Boltzmann equation, assuming a fixed electrostatic potential distribution: 
the potential in turn is obtained by a solution of the Poisson equation assuming a fixed 
charge density distribution. Hence, the solutions of the two equations must be self 
consistent. The solutions of the Boltzmann equation are the phase-space densities of ions 
and electrons, which in turn lead to the respective particle densities through moment 
integrals. When the mean free path is finite, the phase-space density varies along the 
particle trajectories. When the mean free path is infinite, the Boltzmann equation becomes 
the Vlasov equation, where the phase-space density does not vary along a trajectory. Even 
in this case one must still trace the complete trajectory in order to evaluate the phase-space 
density. Hence, the charge density is a complicated functional of the potential distribution; 
that is, it depends on the potential distribution elsewhere in space, and the Poisson- 
Uoltzmann (or Poisson- Vlasov) problem may be said to be nonlocal. 

We are concerned in the present work \ r ith the solution of plasma-sheath Poisson-Vlasov 
systems. In general, one must formulate the problem on a grid and use iteiative numerical 
techniques (Poisson-Vlasov iteration) to obtain self-consistent solutions.* Physicists have 
devised iterative methods for the numerical solution of particular plasma-sheath problems, 
such as those of probes (References 1 to 7), high-velocity satellite wakes (References 8 to 
1 1 ), ion engines (Reference 1 2), and plasma diodes (References 13, 14). Among the 
papers cited, References 5 to 7 and 9 refer to a theoretical analysis of the iteration method; 
these analyses were relatively specialized to the specific problem. Although it is 
perhaps obvious that the various plasma-sheath problems exemplified by References 1 to 
14 are related, their interrelationship and the theoretical basis of the iterative methods 
presented have not been studied systematically. It is felt that such a study would be 
valuable to workers contemplating new sheath calculations, by helping to provide a basis 


♦An alternate approach is that of time-sequential solution (direct computer simulation), although to date this 
has proved feasible only for problems with simple boundary conditions (for example, see papers by R. W. Hockney 
and others in Methods in Computational Physics , 9, B Alder, S. Fernbach, and M. Rotenberg, eds., Academic Press, 
N.Y., 1970). The spherical probe problem has been approached this way by N. Albe r s (Stanford University SU-IPK 
Rpt. No. 499, December 1972) 



for estimating computer costs. Preliminary results of the present investigation have been 
reported in Reference 15. 

It should be noted that mathematicians have investigated the formal theory of iterative 
processes for general elliptic boundary-value problems (References 16 to 21), but have not 
considered their practical applications to the class of problems associated with plasma 
sheaths. 

A unique problem arises when one boundary surface is at infinity, such as in probes 
(References 1 to 7) and satellite wakes (References 8 to 1 1 ). The outer boundary condi- 
tion must be simulated on the computer by a necessarily artificial finite condition (for 
example, an assumed linear condition on the potential and/or its gradient on the surface 
of a large box surrounding the region of interest). This document deals primarily with such 
problems. In it are investigated, both by theoretical analysis and computer experiments, 
two iterative algorithms with respect to stability and number of iterations required. The 
two algorithms are complementary in that the first one depends primarily on the outer 
boundary condition and only slightly on the mesh interval, whereas the second one is only 
siigmiy affected by the outer boundary condition and depends primarily on the mesh 
interval. The first algorithm (Equation (2)) is the prototype of that which has been used 
apparently by most physicists (References 1, 3, 4, 6, 7. 9). but which has only recently 
begun to interest mathematicians (Reference 16). The second algorithm (Equation (7)) 
has received little practical but much mathematical attention (References 18, 19, 22, 23). 

To obtain an adequate quantity of empirical numerical data with a relatively small computer 
expenditure, a one-space-dimensional model problem has been chosen: a spherical probe 
in a collisionless monoenergetic plasma (References 1, 2, 5, 7, 24). Previously, efficient 
numerical methods that take maximum advantage of symmetry (References 1 and 5) have 
been applied to this one-dimensional problem. In this paper the solution is obtained by 
applying a general iterative method that is insensitive to both symmetry and the number 
of dimensions. 

FIRST ALGORITHM 

It is assumed that in an arbitrary three-dimensional problem the partial differential Poisson 
equation has been replaced by a discrete approximation, that is, by a set of difference 
equations based on a chosen grid (and including a chosen outer boundary condition). The 
term solution vector is defined as that vector whose components are the values of the 
potential at the grid points. This solution vector satisfies a matrix equation of the form 

u = f(<T) (i) 

where L is the discrete analog of the Laplacian operator, <j> is the solution vector (potential 
distribution), and F consists partly of the negative of the charge-density vector whose 
components arMhe values of the charge density at the grid points. In^addition, each 
component of F generally depends on more than one component of <p . It is assumed that 
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a numerical procedure is available that is capable of computing F when <j> is given; that is, 
the procedure obtains solutions ofthe Boltzmann equation leading to the charged-particle 
densities and the vector F. Thus, F formally can be considered as resulting from a non- 
linear operation on 0 . * 


The first of tht two iterative algorithms (Equation (2)) will be discussed in detail: 

L0 n+1 =aF(^) + (l^)L0 n (2) 

— ► — ► 

In this algorithm. <p n represents the n-th iterate for the \ector <j>, n denotes the iteration in- 
dex. and a is a positive scalar relaxation parameter between zero and unity that couples (or 
mixes) successive iterates. (More generally, a can be nonstationary (dependent on n), or 
a diagonal matrix (dependent on position) (Reference 1).) When a is 1, there is direct 
iteration, which is the discrete analog of the well known Picard method of successive 
approximations for nonlinear two-point boundary-value problems (References 17, 18, 

20, 2 1 ). When a is small, previous iterates are weighted heavily on the right-hand side of 
(2). It will be shown that the iteration properties of (2) depend strongly on the boundary 
condition. To study the growth of errors, c ■■'ill denote an error vector (defined as the 
difference between 0 n and the true solution vector). Assuming that F is differentiable, 
first-order linearization analysis of (2) leads to a matrix relation between the n-th ana 
(n+ 1 )-th error vectors: 


where the error propagator matrix M is defined by 


M = 1 - a + a L 1 J (4) 

In (4), L * is tne inverse of li'c operator l, and J is the Jacobian matrix of partial 

— ► — > 

derivatives of the components of F with respect to the components of <p . herefore, the 
iteration either converges or diverges depending on whether the spectra! radius of M is 
less than or greater than unity. Because the Jacobian matrix is difficult to represent 
analytically, its effects are studied in a crude manner by introducing an equivalent scalar 
parameter dF/d$. For the situations of interest here, where a must be small, dF/a£ plays 
the role of a Lipschitzian constant (Chapter 9 of Reference 21 or p. 2 of Reference 23). 

This is probably justified as a zero-order approximation if the range of the eigenvalues of 
the Jacobian matrix is small compared with the range of the eigenvalues of L‘* . (Some a 
posteriori support will be given in the section that discusses the second algorithm.) 


* The example of a spherical electrostatic probe is treated in detail in the next section. 
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For simplicity, a three-dimensional Dirichlet boundary-value problem contained within a 
large cube that is divided into N mesh intervals of uniform width h along each of the three 
dimensions is assumed. If the simplest centered seven-point difference scheme is used for 
L, the eigenvalues of L are given hy simple trigonometric functions (for example, by 
generalizing Reference 19 or p. 230 of Reference 22 to three dimensions), and the triply- 
indexed eigenvalues of M become (Appendix A of Reference 6): 


jke 


1 - o - a 




(5) 


(j, k.8-1 N-l) 


According to (5), in the a, X plane, all eigenvalues lie within a fan of straight lines with a 
common vertex (a, X) = (0, 1 ). Let o denote the spectral radius of M, that is. the largest 
eigenvalue magnitude. Then, as shown in Figure 1, a is less than unity and the iteration 
converges for a restricted range of a between zero and a critical value a,, which is controlled 
by the largest eigenvalue of L' 1 . The factor by which the error is reduced after n iterations 
is of the order of 0 ", which for a close to unity is gi.ca approximately by the expression 
(l-n(l-o)]. Fora>a c , a > 1, and the iteration theoretically diverges. The spectral radius 
o(a) i; minimized at the optimum value a , which gives the maximum rate of convergence: 
the minimum number of iterations. In the range between zero and a , the rate of con- 
vergence increases with increasing a. In this range, the dominant eigenvalue is given by the 
largest indices in (5) and is positive; hence, the iteration is monotone (that is, the error 
does not change sign in successive iterations). For a > a , the dominant eigenvalue is 
given by the smallest indices in (5) and is negative; therefore, the iteration is oscillatory 
and the error changes sign in successive iterations. For a > a c , the iteration not only 
diverges but diverges in an oscillatory fashion. 

In problems with only one or two dimensions, the eigenvalues are given by an equation 
similar to (5), but which has only one or two terms in the parentheses instead of three. 

For N » 1 , which is usually the case (Appendix A of Reference 6): 


T dFN 2 h 2 l 
as 2 1+CK — I 

C L 7T 2 J 


( 6 ) 


where C is the reciprocal of the number of dimensions, and K is unity for Dirichlet condi- 
tions. Note that a £ depends on Nh rather than on h (p. 1 59 of Reference 21 ), where Nh 
is a linear dimension across the grid, say in Debye lengths. Hence, if the linear dimensions 
of the grid are large in Debye lengths, so that Nh is large, then (6) shows that the range of 
a within which the iteration converges becomes small, that is, is bounded between zero and 
a c % (Nh)' 2 . In fact, a qualitatively similar relation between convergence and grid dimen- 
sions has been observed by Laframboise (Reference 1 ), Parker (References 6 and 7), and 
Parker and Whipple (Reference 3) on probe problems (Reference 15), and by Call 
(Reference 10) and Maslennikov and Sigov (Reference 9) on the satellite wake problem. 
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Figure 1. Spectral radius of iteration matrix for linearized problem. 

Also, if the charge density has a singular character, as in the spherical probe problem to be 
discussed later, dF/d <p can become large, thus forcing the reduction of a still further. 

Finally, it is observed that as (dF/d0)(Nh) 2 becomes large, the following conditions prevail: 

• a . tends to coincide with a while both become small. 

opt c 

• The iteration is monotone, ard, for a within the convergent range, the spectral 
radius of M is controlled by the smallest eigenvalue of L' 1 and is given by 1-Aa 
where A = 1 + C (dF/d0)(h 2 /4). 

• Because Aa « 1, the reduction of error after n iterations is approximated by 
the expression exp ( -Ana). 

• For a even slightly larger than a c , the iteration diverges rapidly. 

Incidentally, Lettenmeyer’s criterion (References 17, 18, and pp. 40 and 143 of Reference 
21 ; the authors were not able to locate Lettenmeyer’s original paper) for the convergence 
of the Picard procedure in the one-dimensional continuous Dirichlet problem follows 
immediatelv from (5) with a = 1 and N » 1. That is, the spectral radius is then given by 
N 2 h 2 |dF/d0|/jr 2 , which must be less than unity for convergence. This is the Lettenmeyer 
criterion, where |dF/d0| is actually a Lipschitz constant. 

To evaluate the role played by boundary conditions, a two-surface boundary-value problem 
with a Dirichlet condition on one boundary surface and a floating linear relation between 
the potential and its gradient on the other boundary is considered. This has been employed 
in probe problems to simulate the condition at infinity (References 1, 3 to 7, 1 5). The 
largest eigenvalue of the reciprocal matrix L' 1 is increased by a factor K > 1 which appears 
in (6). '♦ can be shown, by analyzing the roots of the appropriate determinantal equation 
(for ex, ,.iple, Equation (1.34*) of Reference 25) for a one-dimensional problem with a 
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Dirichlet condition (0 = 0) at one end and a floating condition ( <p ' ~ -k$) at the other end, 
that K = 7r 2 /z 2 , where z is the root of kNh = ztan(z-7r/2) lying in the range ir/2 < z < 7r. 

That is, K goes from unity to 4 as k goes from infinity (0 = 0) to zero (</>' = 0). Thus, the 
use of a floating condition tends to require smaller values of a and correspondingly more 
iterations as opposed to the Dirichlet condition. The value of dF/d0 may change in such a 
direction as to cancel this effect. It has been observed in numerical experiments with the 
spherical probe (Table 1 ) that in some cases floating conditions converge more rapidly, 
and in other cases more slowly than fixed Dirichlet conu.tions. The value of dF/'d<£ is not 
usually predictable and is considered to be an empirica’ly determined parameter. Never- 
theless, it is a useful parameter because, from problem to problem, it varies more slowly 
than a c (Table 2). 

The following numerical experiments with the spherical probe model also sne aid be rele- 
vant to three-dimensional problems. Preliminary results were reported in Reference 1 5. 

SPHERICAL PROBE 

The problem of the spherical electrostatic probe is used as the model for computational 
experiments involving tests of three types of numerical boundary conditions simulating 
the condition at infinity, and two types of Poisson-Vlasov iteration (the first and second 
algorithms). In particular, the monoenergetic version of the probe calculation is used, 
which assumes that all attracted particles have the same nonvanishing total energy E 0 . 
Theoretical calculations of current-voltage characteristics for this model, with repelled 
particles assumed to be described by a Boltzmann factor (that is. an exponential function 
of local potential), have been made by Bohm et al. (Reference 24), Bernstein and 
Rabinowitz (Reference 2), and Laframboise (Reference 1). A similar model, but one 
using repelled as well as attracted monoenergetic particles, has been used by Guderley 
(Reference 5) (based on unpublished work by G. Medicus), and is the model discussed in 
this section. These monoenergetic models are physically justified on the basis that the 
problem can be solved numerically for arbitrary energy distributions (usually a Maxwellian) 
by summing over monoenergetic contributions (that is, by energy-quadratures). Such 
calculations have been done independently by Laframboise and Parker (unpublished work) 
duplicating results obtained earlier by Laframboise (Reference 1), who used a special method 
for Maxwellian energy distributions. Another justification is based on computed solutions 
that show that as the probe potential becomes large the monoenergetic attracted-particle 
current becomes identical to the Maxwellian current. 

Figure 2 indicates the parameters of the probe problem, with <f> denoting dimensionless 
potential in units of E 0 /e, where E 0 is the singular energy (the same for both attracted and 
repelled particles, consisting of electrons and singly-charged ions with charge e). The unit 
of length is (E 0 /2ffn 0 e 2 ) l/2 , as in Reference 5, that is, larger than the usual Debye length 
by a factor of y/T. The range of r is between r p (dimensionless probe radius, where the 
potential is 0 p ) and r b , the dimensionless outer-boundary radius. 


6 




Figure 2. Spherical collisionless monoenergetic probe model: Artificial 
boundary conditions are established to simulate boundary condition at 
infinity (r^ is moved outward until current does not change). 

The equations that may be used for a monoenergetic problem in which the ions and elec- 
trons both have the same energy are presented in condensed form as follows. (See 
References 5, 7, and 24 for derivations and more detailed discussions.) 

In the Poisson equation 

d 2 <t> 2 d <t> 

= p(0) = 2n (-0) - 2n(0) 

dr r dr 

where n(-0) and n(0) denote dimensionless electron and ion density, respectively. These 
densities are functionals of a function g(r), defined by 

g( r ;= r 2 (l-0) 


where g(r) is proportional to the square of the largest angular momentum or impact 
parameter which a particle can have at p^nt r, and is such that there is a turning-point at 
r (Section 4.2 of Reference 24). Minima in g, namely gj and g m , correspond to maxima 
in the effective potential and govern the density and current as seen in the following 
equations: 



g = Min gfr) = Min g, (r) 

r r 
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• If g(r) > 0 for r p < r Q < r < r b , then 
^=0 

n(0) = 0 for r < r Q 



• Dimensionless current density is 

j B 

— = (j = current at zero probe potential) 

Jo r P 

Finding the minima and g m (where g x is the lowest value of g ahead of a Fixed point r, 
and g m is the minimum over all r) is an operation easily performed and ideal for computer 
calculations when the 0-values are defined at a set of grid points. The special radius r Q is 
defined as the smallest radius such that g is positive for all r > r Q , and is the largest radius 
for which 1 - 0 assumes negative values. In ai! of the cases treated here, the g{r) for the 
correct solution has a single minimum at the absorption radius r c as shown for the 
attracted particles in Figui 3. A typical function g(r) for the repelled particles is also 
shown in Figure 3 by the dashed curve, where there is a radius r Q such that g is negative for 
r < r Q . For intermediate potential distributions occurring during iterative cycling, the 
function g(r) may become very complicated, but the above formulas take into account any 
possible variation. The formulas fo* monoenergetic particles are directly useful for arbi- 
trary polyenergetic distributions by means of simple quadratures. 
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Figure 3. Special case of a nonmonotonic function g(r) with a 
single minimum at r c for attracted particles. For repelled 
particles, g(r) may change sign (as at r Q in dashed curve). 




The quantity of interest is the current of ions or electrons collected by the probe. Since 
r b cannot be infinite, the question is how large must r fc be in order to simulate sufficiently 
weil the boundary condition at infinity where 0 vanishes. It is found in practice that as r b 
increases, the current wil! eventually level off and become independent oi r b - It is also 
found tnat increasing the value of r b causes the computer cost per solution to incrc >se, 
since mare grid points would be required to control the truncation error, and the number 
of iterations would tend to increase in accord with the analysis of the preceding section. 

Several different types of boundary conditions may be imposed at r fe . ( i ay lor (Reference 
26) has studied various types of conditions with respect tc how well they simulate the 
boundary condition at infinity; the effect of the condition on iterations was not studied.) 

Three types of boundary condition studied in the present work are illustrated in Figure 2; 
(a) the zero-potential fixed condition, which has been used aa hoc for satellite wake calcu- 
lations (References 8 to 11); (b) tne inverse-squareuaw floating condition, which represents 
the proper asymptotic behavior lor the probe problem, and so has been used effectively 
for spherical probe calculations (References 1 and 5), and (c) the zero-gradient floating 
condition. The zero-potential and zero-gradient conditions represent opposite-extreme 
conditions. General liuear relations, including inverse-sq M are and exponential laws, should 
give results lying between these extremes. Therefore, the three conditions studied are 
probably representative of the range of interest, t hese conditions have different properties 
with respect to the leveling-off value of r b and the number of iterations required for 
convergence; it is of interest to determine which condition results in the smallest net 
computer time per solution Therefore, the optimum values of r b and the iteration 
parameters for both the first and second algorithms will be determined. 

Four physical problems have been studied: probe radius r p equal to 10 and 100 uni 
(14.1 and 141 Debye lengths), with probe potential 0 p equal to 10 and 100 units f m 
probe radius. The extreme values of probe radius and potential were chosen because they 
are physically interesting and also numerically interesting, that is, likely to be costly in 
computer time. For each problem, all tliree boundary conditic s were investigated to find 
the smallest boundary radius r b consistent with stationary current, and the smallest number 
of iterations required for the first algorithm (Equation (2)). To compare numbers of 
iterations, the iteration was always started with the Laplace solution as a standard initial 
estimate. The iteration was continued until successive iterates, if they converged, differed 
by less than one part in 10 s in the current density 

Resul s for the 12 cases ate shown in Table 1. For each case, the table gives the values of 
four quantities: the minimum r b for stationary current, the optimum (essentially the 
critical) value of a, the minimum number of iterations n min , and the approximate time in 
seconds on an IBM 7094 computer. (Empirically, the time was about one ms/grid point/ 
iteration.) For all but one case the mesh interval was b=l, yielding the values shown in 
Table 1 of dimensi nless current density j/j Q , with a truncation error of the order of a 
few tenths oi one percent, as indicated by the numbers in parentheses. For the zero- 
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potential boundary condition and for r p = 0 p ~ a the smallest value of h used was 2.5 
because of computer-time limitations. The truncation error was estimated by running 
with successively smaller values of h (down to 0. 1 ). The computer time in each case is 
proportional to the product of n and the number of grid points. (This is because a direct 
elimination method was used to solve the tridiagonal system of Poisson difference equa- 
tions, so that the time per grid point is a constant for that part of the calculation.) 

The boundary conditions may be assessed with respect to computer economy as follows. 

• For 0 p * 10 (large but not extremely large potential), the two floating conditions 
are comparable. 

• For = 1 00 (extremely large potential), the inverse-square floating condition is 
superior to the zero-gradient floating condition. 

• For r p = 10 (large but not extremely large radius), *he fixed condition is 
comparable or superior to both floating conditions. 

• For r p = 100 (extremely large radius), the floating conditions are superior 
to the fixed condition. 

For large probe potential or probe radius, the fixed condition has the peculiar property 
that, for a fixed a, as r b increases, the number of iterations goes up and then down again. 
For Case 10, as noted in Table 1 , convergence was essentially impossible to obtain when 
r b was in a certain range of values. The floating conditions were free of such difficulties. 

Table 2 presents the empirical values of n/(Nh) 2 , a(Nh) 2 , K, and (dF/d0) obtained fi^m 
(6), where Nh = r b - r . Throughout the table, the ratios of the largest to the smallest 
values of n/(Nh) 2 , a(Nh) 2 , and (dF/d<£) are 9, 15, and 12, respectively. These factors are 
much smaller than the variation of a in Table 1, that is to say, by a factor of 300. For 
9 out of the 1 2 cases, dF/d 0 has values between about 1 and 3. In Cases 3, 6, and 10. 
dF/d0 has values of about 5, 7, and 16, respectively. The occurrence of these larger 
values is probably connected with the existence of a singularity in the derivative of the 
density as a function of potential (where a radicand vanishes). 

If the error reduction by iteration is described by the expression exp ( -Bna), the em irical 
coefficient B is found sometimes to be less than unity, n> contradiction to the linearization 
theory. This is probably due to the closeness of the spectral radius o to unity, so that 1-a 
is sensitive to ~~all perturbations in the theory. For example, the neglect in the spread of 
the eigenvalues of the Jacobian matrix J could lead to such a perturbation. This affects 
the number of iterations. However, the critical value of a is much less affected. Indeed, 
consistent with the theory, a, is found to be a decreasing function of Nh, and is insensi- 
tive to h for a given Nh. 

Aside from the properties of the iteration process, Table 1 shows that the minimum value 
of r b is larger for the fixed condition than for the two floating conditions, but is independ- 
ent of which floating condition i used. This is consistent with References 1 and 26. 
Moreover, the minimum r b is found f o be independent of h (not indicated in Table 1 ). 
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As a result of the above findings, the toilowing procedure is recommended in order to save 
computer time when solving similar problems for any chosen boundary condition: First, 
minimum values of r b and critical values of a should be found while using large values of 
h. that is, few grid points. Then h can subsequently be reduced in order to diminish 
truncation error without the necessity of further changes in r h and a. It is not yet clear 
how to choose the proper boundary condition to be used for an arbitrary problem. 

An interesting but unjustified assertion has been made by Maslennikov and Sigov 
(Reference 9) in connection with their use of (2) to iterate the high-velocity satellite 
wake problem. To ensure stability, the volume of the grid (the region cf interest) must be 
smaller than a number proportional to the ion energy. Tne present analysis defined Nh 
to be an appropriate linear dimension of the grid in Debye lengths, and found that (Nh) 2 
must be iess than a certain constant for a fixed a. If X represents the appropriate linear 
dimension of the grid, and if the effective Debye length is based on the ion directed 
kinetic energy E Q (rather than on the electron thermal energy kT). then (Nh) 2 is propor- 
tional to X 2 /E q . Hence, the criterion is that a cross-sectional area (X 2 rather than the 
volume X 3 as stated in Reference 9) of the grid must be smaller than a number propor- 
tional to E q : The characteristic energy is the larger of E Q ana kT. Thus, the criterion for 
stability requires that the dimension X be smaller than a certain number of effective 
Debye lengths A*, where A* is defined by \/E /4irn 0 e 2 when E 0 » kT (Reference 9), 
and by\/kT7<irn 0 e 5 when E 0 $ kT. 

SECOND ALGORITHM 

This section discusses the second algorithm: 

(L-p)« n+1 =F(0 n ) - p? n (7) 

where p is a scalar relaxation parameter. This extension of the Richardson method (p. 226 
of keterence 11 , p. 90 ot Reference 23) has been studied in connection with mildly 
nonlinear elliptic equations (References 1 8 and 19). It will be seen that the iteration 
properties of (7) are essentially independent of the boundary condition, in contrast to (2). 
This is because the eigenvalues of L, which contain the boundary condition information, 
have a minor effect on the spectral radius of the iteration matrix M that relates successive 
error vectors. By applying a linearization process to (7), similar to that applied to (2), 

M = (p-L)-' (p-J) (8) 

where J is the Jacobian derivative matrix as in ''4). (If the scalar parameter p in (7) were 
replaced by the matrix J, the lin ar iteration matrix M would vanish, resulting in a Newto- 
nian-type process, that is, a disc’^te analog of quasi-linearization (Reference 20).) It is clear 
that the eigenvalues of J play a crucial role in determining the eigenvalues of M. Assuming 
the effect of L to be that of a perturbation, then to lowest order the eigenvalues of M are 
approximated by 



( 9 ) 


X - — 

P -x L 

where Xj is an eigenvalue of the Jacobian matrix J, and X L is an eigenvalue of the operator 
L. 


It is assumed that the eigenvalues of J are all positive* and lie in the range (X t < \ } < X 2 ). 
Also, under the conditions leading to (5), all of the X L are negative. Hence, for fixed p and 
Xj, the largest value of |X| occurs for the minimum value of -X L , which is given by 
X 3 = *r 2 /(KN 2 h 2 ) when N » 1. Thus, if p is large compared with X 3 , the spectral radius 
of M, namely a, becomes independent of the eigenvalues of L and is approximated by 


G 



for p < 


X 2 + X, 


( 10 ) 


forp>^^ ("> 

Note that, considered as a function of 1/p, o(l/p) is qualitatively similar to o(a) plotted 
in Figure 1 . The minimum value of o, corresponding to the minimum number of 
iterations, is 


a 

m,n X 2 + X, 

which occurs when p is equal to the optimum value 


( 12 ) 


n = 
A opt 


X 2 + Xi 


(13) 


When p exceeds the critical value p c the spectral radius a exceeds unity, and the iter on 

rl m ■ rv»rri,t' nrViurp 


D 
• C 



(14) 


Therefore, if the eigenvalues of J are approximately equal (that is, X, ~ X 2 ), then 
o . « 1 and p ~ 2p ~ X,. However, if X, « X,, 

min r opt r c J 1 2 



(15) 


Moreover, if o is close to unity, then according to (1 1 ) the error after n iterations is 
reduced approximately by the expression exp ( -nX! Ip). 


•Just ' table on the basis of analytic forms approximating the space-charge density (Appendix A of Reference 6). 
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Thus, empirical evaluation of p opt and p c from iteration experiments with our model 
problem will give information regarding the eigenvalues of the matrix J for the s roblem. 

A check can also be made on the assumptions leading to (5) and (9), and a direct 
comparison of the two algorithms can be made. Presented in Table 3 are the results of 
computations using the second algorithm for the same 12 problems presented in Table 1 
for the first algorithm. The minimum values of r fe are the same as in Table 1 and are not 
shown in Table 3. For each case, the table presents five quantities: the optimum value 
p t , the range of p containing the critical value p^. the minimum number n of iterations, 
the approximate value of X 2 , and the value of X 4 = 1/X 3 (the dominant eigenvalue of -L' 1 ). 

It is evident from Table 3 that p t , p c , and the minimum number of iterations are inde- 
pendent of the boundary condition. Actually, for a given r p and 0 p , they depend only on 
the mesh interval h. Hence, for a given truncation error, governed by h, the two floating 
conditions require equal numbers of grid points but fewer points than the fixed condition; 
therefore, the floating conditions are equally economical and are more economical than 
the fixed condition. 

In all cases, the fact that the fractional difference between p opt and p c is very small is 
evidence of the validity of (15) and of the inequality Xj « X 2 (assuming (9) to (1 1) are 
valid). Hence, from (15) and the empirical values of p , a 2 may be estimated. Now 
P 3pt » X 3 , where X 3 = 1/X 4 , in all cases. Therefore, the assumption leading to (9) to 
(11) is justified, that is, that the effect of L is that of a perturbation. Also, the fact that 
X 4 is in all cases larger than X 2 supports the assumption that the range of the eigenvalues 
of J is small compared with the range of the eigenvalues of L 1 These considerations 
apply to h = 1. However, as noted in the footnotes of Table 3, when h was set equal to the 
small value 0. 1 instead of to 1 , it was found in Case 1 that p . was reduced to 0.02, a 
considerable reduction. Because p op< « X 3 in this case, the effect of L can no longer be 
a small perturbation, so that (9) to ( 11 ) are probably invalid. However, the number of 
iterations was 47, which is approximately the same as for h = 1 , which was 52. 

C' iparison of the numbers of iterations in Table 3 with the corresponding numbers in 
.able 1 shows that for the given values of h, the second algorithm is more economical 
than the first. In some cases, the number of iterations is even reduced by an order of 
magnitude (for example. Case 10). However, the second algorithm is so sensitive to the 
value of h that the situation car. be significantly reversed. For example, for Cases 10 to P 
in Table 3, no satisfactory convergence was obtainable for h = 2.5; that is, the behavior 
was peculiar in that for p = 90 the current appeared to converge after 5000 to 6000 itera- 
ti ms, but to a value 4 percent lower than that resulting from the first algorithm. For p on 
eit* ;/ side of 90, however, the iteration number was greater than 25 000, the maximum 
allowed. Nevertheless, in all cases shown in Tables 1 and 3, the iteration was assumed to 
have converged properly because both algorithms converged to the same value in every case. 

Goddard Space Flight Center 

National Aeronautics and Space Administration 
Greenbelt, Maryland February 19, 1973 
188-48-52-05*51 
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